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ABSTRACT 


The traditional study of underwater explosions (UNDEX) with respect to ship 
damage became of interest during World War II when torpedo explosions near a ship 
created more damage than a direct hit. Following the war, many full scale ship shock 
trials were conducted that provided much of the empirical data that is used in the field 
today. However, one type of shock phenomena became of interest in the late 1960s that 
potentially could be more damaging than a typical underwater explosion; an implosion. 
Crude implosion experiments were conducted in the late 1960s. Although these 
experiments collected data on pressure waves, more emphasis was placed on the 
acoustical properties associated with an implosion event. Today, one of the Navy’s 
concerns is about the potential for the implosion of a pressure vessel in close proximity to 
a submarine hull. A computational approach is desired that will predict the source 
strength of an implosion. 

This thesis will cover the basic principals of underwater shock phenomena, 
including explosions and implosions. Drawing from previous experiments and 
computational simulations, a detailed investigation of the implosion event will be made 
using, DYSMAS, a coupled Eulerian-Lagrangian solver. DYSMAS will be used to 


compare the characteristics of implosion and explosion events. 
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I. INTRODUCTION 


A. BACKGROUND 

The field of underwater explosion (UNDEX) research is not a new concept. 
During World War II, it was noted that an UNDEX event could actually damage a vessel 
without directly impacting the target with a torpedo or mine. The concept of “close is 
close enough” was especially applicable to submarine warfare. In particular, depth 
charges were used to destroy submarines during this era. A depth charge did not 


necessarily need to physically touch the submarine in order to cause catastrophic damage. 


Numerous studies were conducted on both surface ships and submarines after the 
war. Many of these studies focused on the vulnerabilities of ships and submarines to 
UNDEX phenomena. At the time, the studies focused on the initial pressure wave 
created by the explosion. The subsequent pressure pulses caused by bubble oscillation 


were also studied, but to a lesser degree. 


Although less emphasis was placed on the study of the bubble oscillation, one fact 
remained apparent. The initial shock wave indeed excited the response of the ship or 
submarine, but intuitively, it does little for the low frequency response of the ship. 
Something else was causing the low frequency vibration of the hull girder. The 
explanation for this low frequency vibration, or whipping as it would be later dubbed, 
was indeed the oscillating bubble underneath the hull of the vessel. Researchers quickly 
realized the effect of whipping, which could possibly “break the back” of the ship or 


submarine [Ref. 1]. 


The cause of whipping is not necessarily the bubble oscillating, but the 
subsequent pressure gradients cause by the bubble oscillation. As previously stated, the 
initial shock wave provides the energy needed to excite many of the natural frequencies 
of the ship or the shipboard equipment. However, it is the gas bubble oscillation that 
remains the underlying cause of low frequency excitation [Ref. 1]. Whipping occurs to 
the extent that it does for two main reasons: 1) the bubble oscillations occur at a 


frequency that is roughly equal to the low mode natural frequency of the hull girder, thus 


creating a resonance condition, and 2) because of the long duration flow loading from the 
initial bubble expansion which initiates the ship movement. 
B. SCOPE OF RESEARCH 

It has been established through testing that explosions near a submarine or surface 
vessel create a damaging initial pressure wave and a subsequent equally as dangerous 
oscillating gas bubble. However, an explosion is not the only method of producing a 
pressure wave and oscillating bubble. Perhaps an implosion event near a submarine 


would also create a damaging pressure wave and bubble phenomena. 


For instance, suppose a submerged submarine was a host vessel to a smaller 
experimental pod or pressurized delivery vehicle. The interior compartment of the 
vehicle would be at atmospheric pressure, which would be much less than the pressure of 
the surrounding hydrostatic environment. If the vehicle were to encounter damage in 
some way or perhaps even randomly implode in close proximity to the host submarine, it 
would indeed create a pressure wave and oscillating bubble action that would ultimately 


impact the host vessel. 


Computer modeling of implosion events are now becoming of interest. During 
the late 1960’s, experiments were conducted on the implosions in pressure vessels. 
These experiments indicated that the maximum hydrodynamic pressure decreased as the 
distance from the implosion center increased [Ref. 2]. In other words, if an implosion 
event occurred in close proximity to a submarine, the vessel would encounter a large 
pressure wave. This pressure wave is a result of the surrounding water rushing inward 
upon the collapsing structure and low pressure air inside the structure. The surrounding 
water builds momentum as it rushed inward during the collapse. When the air reaches a 
minimum volume, the velocity of the water is forced to zero and the water compresses, 


resulting in a shock wave that travels back out into the water. 


This thesis will investigate the pressure wave and gas bubble oscillation caused by 
an implosion event. The research conducted will utilize an Eulerian-Lagrangian coupled 
solver known as DYSMAS. In addition to studying the pressure wave cause by an 
implosion event, a comparative study will also be done to identify similarities between 


explosion and implosion events. 


I. UNDERWATER EXPLOSION THEORY AND MODELING 


A. ORDER OF EVENTS FOR AN UNDERWATER EXPLOSION 

An underwater explosion releases a compressive shock wave which begins to 
move at 25,000 ft/sec and exponentially slows down to nearly 5,000 ft/sec (which is the 
speed of sound in water) as it moves outward. In order to simplify the mathematics of the 
problem, it can be assumed that the shock wave travels at the speed of sound in water, 
5,000 ft/sec. By using this approximation, the following formula for a plane wave 


approximation in a compressible flow environment is: 


P=p,C,u 


where P is the fluid pressure, Co is the speed of sound in water, po is the density of the 


fluid, and wu is the fluid particle velocity. 


Over time, empirical formulas have been developed to describe the fluid particle 
velocity and pressure during the initial shock phase of an UNDEX event cause by a 
specified charge type and weight. The following equation is an example of one of the 
empirical equations formulated to determine peak pressure at a specified distance from 


the charge: 


w’) 





Pe, = K, : 
in which case Pingx is the peak magnitude of the pressure in the shock front (measured in 
psi), K; and A; are parameters which are dependent on the type of charge (TNT, HBX-1, 


etc.), and R is the radial distance from the explosive charge to the target (measured in 


feet) [Ref. 3]. 


As the shock wave propagates through the water, it will inevitably reach either its 
target (a ship or submarine) or reflect off of the seafloor or the air-water interface. Shock 
waves that reflect off of the seafloor primarily remain compressive and reflect or transmit 


in accordance with Snell’s Law, which states: 


C C 
sina sina’ sina" 


0 0 0 





where Co is the speed of sound in a given environment while the three angles are shown 
in Figure 1 below. 
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Reflected Shock Front Velocity - Incident Shock Front 
Figure 1. Illustration of Snell's Law (From Ref. 3) 


The shock waves that reflect off of the air-water interface are known as surface 
reflection or rarefaction waves. These waves are tensile in nature. In addition, these 
reflection waves also follow Snell’s Law. The direction and magnitude of these waves 
are determined by imagining an “image” charge which is equidistant of the charge below 
the waterline. Figure 2 illustrates the “image” charge and Snell’s Law for both bottom 


and air-water interface reflection waves. 
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Figure 2. Typical UNDEX Phenomena Geometry (From Ref. 3) 


Notice that the tensile wave as a result of surface reflection creates a tension force 
in the water. This tension causes a foamy appearance at the surface of the water. This 
phenomenon is known as bulk cavitation, where the pressure of the water is essentially 


equal to zero. Figure 3 is a photograph of this phenomenon. 





Figure 3. Photograph of Bulk Cavitation Due to UNDEX Event 


After the initial shock wave passes, subsequent pulses are experienced at time 
intervals that are functions of the type of explosive, charge weight, and depth of the 
charge. When an explosion occurs gasses are released. These gasses form a large 
expanding bubble which continues to expand until hydrostatic pressure halts the bubble’s 
growth. When the bubble is finished expanding, the pressure inside the bubble is much 
less than that of the outside hydrostatic pressure. This causes the bubble to quickly 
collapse. When this bubble collapses, it releases a “bubble pulse”. This pulse in turn 
creates another expanding bubble of lesser magnitude. The term “shock wave” is 
replaced with “pulse” because the pressure waves are now being created by the 
oscillating bubble. Recall that the initial shock wave, which causes initial damage to a 
ship or submarine, is of short duration and high magnitude. However, the pressure waves 


caused by the bubble pulses are of lower magnitude and longer duration. 


Due to hydrostatics, the oscillating gas bubble migrates to the surface at varying 
speeds. As the bubble size decreases, the bubble moves more quickly to the surface. On 
the contrary, as the bubble re-expands, its velocity decreases. This varying velocity is 
due to the drag forces on the bubble. After many bubble oscillations, the energy from the 
explosion is either dissipated into a large number of small gas bubbles or the bubble 
reaches the surface and vents to the atmosphere [Ref. 4]. As the bubble moves toward 
the surface, the bubble dynamics can cause different disruptions to the water’s surface 
depending on the initial depth of the charge. If the explosion is not at great depth, a large 
plume of water may occur at the surface of the water. 

B. MODELING OF THE GAS BUBBLE OSCILLATIONS 

As previously discussed, the initial damage done to a ship or submarine is due to 
the initial shock wave. The following pulses which are due to the oscillating bubble 
create the low frequency excitation. This low frequency excitation is known as 
whipping. In order to better understand whipping, one must first understand these gas 
bubble oscillations. Once theses oscillations are understood, mathematical models can be 


created and employed to predict the damage done by the bubble pulses. 


Each bubble collapse results in a pulse that can be modeled as an acoustic wave. 
However, the pulses are not the primary means by which whipping occurs. The low 
frequency mode of the hull girder is primarily excited by the inertial movement of water 
that occurs when the gas bubble oscillates. Thus, whipping is not directly caused by the 
pressure waves associated with the pulses, but by the “fluid flow associated with the 


entire pressure field” which was in fact caused by the gas bubble oscillations [Ref. 5]. 


The motion of the gas bubble oscillations and the resulting pressure effects are 
seen in Figure 4. The vertical positions of the bubble are displayed with respect to time 


which is on the horizontal axis. Note that the bubble does not migrate to the right. 
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Figure 4. Bubble Oscillation and Subsequent Pressure Effects (From Ref. 3) 


There are very important differences between the initial pressure wave and the 
subsequent gas bubble pulses. The distinctions between the two resulted in two 
approaches to mathematical modeling. As previously stated, the initial shock wave is 
modeled as an acoustic wave, which required the assumption of a compressible fluid 
medium. On the other hand, the bubble oscillations are modeled in an incompressible, 
irrotational, and inviscid fluid medium. The physical processes which occur for the 
initial explosion and the smaller “explosions” associated with bubble collapse are 
basically the same, only differing in magnitude and duration. The two mathematical 
models only partially capture what is occurring and are used based on the most dominant 
effects upon the ship at the time frame that is to be analyzed. The model for the gas 
bubble oscillations must capture both the pushing and pulling on the ship caused by the 
pressure gradients due to bubble oscillation. Therefore, a study of the field immediately 
surrounding the gas bubble is needed. Fortunately, the assumptions of irrotational, 
incompressible, and inviscid flow allows for straightforward calculations. However, the 
model of the initial shock wave must capture the large pressure impulse experienced by 
the ship. For this situation, the compressible fluid, acoustic model works well. An 


incompressible fluid dynamics algorithm could be employed to model the initial shock 


wave, but the assumption of incompressibility would be invalid and would result in a 


failure to model the shock wave’s large amplitude. 


Chertock’s original algorithms which modeled gas bubble oscillations made use 
of the assumption of an incompressible, irrotational, and inviscid fluid medium [Ref. 6]. 
When Hicks developed his modeling approach, he utilized the same assumptions. Hicks 
stated that the incompressibility assumption was valid because “the pressure in the gas 
bubble is very low during most of its pulsation cycle...” [Ref. 7]. The pressure is high 
only during the short period of the initial shock wave and also, but to a lesser extent, 
when the bubble collapses. 

1. The Hicks Bubble Model 

The derivation to follow on the gas bubble behavior from an UNDEX event was 
used by Hicks in his models, but was based on derivations by Herring and Taylor in 
addition to the work done by R. Cole which was completed shortly after World War II 
[Ref. 5, 7]. 


The premise of Hicks’ model is the assumption of a fluid point source of strength, 
denoted as “e”, which originates at the center of the explosion or bubble. This source 
approximation is used to model the gas bubble as an expanding sphere. This source 
refers to a volume flux, denoted as “Q’”, coming from a point source. The flow will 
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spread out radially and at a specified radius, denoted as “r’, will be: 
e 
= = 2 2 4 
4mr° or 
u, =0 





where u, and ug are radial and tangential particle velocities, respectively. Thus, it follows 


that: 


The velocity potential is derived as: 


In order to simplify the mathematical modeling, Hicks used the potential flow 
theory for a two dimensional case, using an infinitely long cylinder as the target. Using 
this approximation, Hicks then assumed that the fluid particle velocity surrounding the 
cylinder was uniform about the cylinder’s cross section. This approach was entitled the 
“distant flow” method because it is based on the assumption that the uniform flow 
occurred at a distance “several ship diameters away from the charge” [Ref. 7]. Further 
studies by Hicks revealed that his two dimensional algorithms were indeed accurate for 
explosions that were “as close as one-fifth of the width of the hull transverse to the 


direction of vibration” [Ref. 8]. 


Equations 4 and 5 represent expressions that model a stationary oscillating 
bubble. In order to include bubble migration, or the bubble’s tendency to rise to the 
surface due to buoyant forces, Hicks added the dipole term to the velocity potential 


equation, thus becoming: 


é é 
g=—+—-cos(G,) 7 
Bee 
where the first term is the simple source term and the second term is the dipole term [Ref. 
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When Hicks initially allowed for bubble migration, he noticed that his algorithm 
overpredicted the bubble migration rate. The root of the problem was that Hicks modeled 
the bubble as a sphere. For a deep explosion, this assumption is quite accurate. 
However, when an explosion occurs at a more shallow depth, the bubble takes the form 
of a mushroom cloud, which is denoted in Figure 5. The mushroom cloud shaped bubble 
rises more slowly than a spherical shaped bubble. In order to correct this problem, Hicks 


incorporated a high drag coefficient of 2.25 [Ref. 5]. This drag coefficient is larger than 
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those typically associated with a submerged spherical object. However, it has calibrated 
for shallow weapon-sized charges which were based on test data from underwater 


explosions. 





Figure 5. Photograph of Bubble Pulse (From Ref. 3) 


In addition, the Hicks Bubble Model does not account for the dissipating 
amplitude of the bubble pulses. Therefore, the model is only accurate for the first period 
of bubble oscillations. Hence, the excitation of the hull girder due to bubble oscillations 
is stopped after 1.5 bubble oscillations, at which time the girder undergoes free vibration 
[Ref. 9]. Initially, Hicks believed that his modeling technique was inadequate due to the 
aforementioned approximations. However, after further research, Hicks determined that 
the effects on a ship due to gas bubble oscillations were well represented by his initial 
model. He found that his distant flow approximation actually proved to be quite adequate 


at normal standoff distances. However, when a charge is placed near a target his 


11 


approximation loses merit. In Hick’s defense, most other approximations also prove 
inadequate at close range. At a very small standoff distance, the bubble tends to engulf 
the hull which causes localized damage. This damage will also tend to further distort the 
flow around the target [Ref. 1]. 
c. THE PHYSICS OF COLLAPSING WAVES 

The implosion phenomenon is also an important part of UNDEX theory. In the 
past, much work has been completed concerning explosion produced shock waves and 
subsequent pressure pulses caused by bubble oscillations. However, significant increases 
in the initial pressure wave can be accomplished by using an implosion creating device. 
Because the damage producing ability of an underwater detonation is more influenced by 
the initial shock wave rather than the pressure pulses caused by the bubble oscillation, it 
is desirable to increase the amount of energy that contributes to the initial shock wave. 
Of course, this is done at the expense of smaller pressure pulses. One method for 
producing such an increased initial pressure wave is shown in Figure 6. Notice that the 
detonation waves are converging on the center of the sphere, rather than diverging as 


done in a typical underwater explosion [Ref. 10]. 
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Figure 6. Solid Explosive Charge Surrounded with Plural Detonators (From Ref. 
10) 
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In Figure 6, a solid explosive charge is provided with many detonators attached to 
the outside surface. At a specified time, divergent waves emerge from the detonators and 
converge toward the center of the charge. The shock waves form an uninterrupted wave 
which closes in on a decreasing volume. Thus, the result is a peak implosion pressure 
within a body of water. Figure 7 displays how the implosion shock wave can be 


significantly greater than a shock wave induced from a typical underwater explosion 
[Ref. 10]. 
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1. The Collapse of a Bubble in a Liquid 

In order to understand the physics behind an implosion, it is useful to consider the 
collapse of an empty spherical cavity in a liquid. Bubbles can form if the pressure at a 
specific point in a flow field is less than the vapor pressure of the liquid. This 
phenomenon is especially prevalent when dealing with cavitating propeller blades of a 
ship or submarine. When these bubbles move at a later time to a region of greater 
pressure, they collapse. The collapse of these bubbles result in a strong shock wave that 


causes damage to surfaces, such as the aforementioned propeller blades. 


In studying the collapse of a bubble in a liquid, one particular approximation must 
be made. This approximation deals with the matter inside of the collapsing bubble. In 
reality, the inside of the bubble is not empty; it contains a vapor. For the calculation of a 
collapsing bubble, the vapor inside of the bubble is neglected. This assumption is 
especially valid in the final stages of collapse when the pressure build up near the cavity 
wall is high. In addition, it is also determined that the bubbles complete their collapse 


when r = 0 [Ref. 11]. 


At an initial time of t = 0, a spherical cavity of radius Ro is in an infinite expanse 
of fluid at a surrounding pressure of Po. The pressure inside the cavity is zero. Thus, the 
pressure difference is equal to Po. Because of the pressure difference, the bubble begins 
to shrink. For the initial motion, compressibility effects are neglected because the motion 
is slow and the pressure rise is low. Hence, the continuity equation for an incompressible 


fluid yields: 


Vinpe + =0 8 
Or r 


because the flow is spherically symmetric. Solving Equation 8 with the boundary 


condition that the particle velocity, u, at the cavity wall, R(t), is equal to the velocity of 
the interface, R (t), the velocity distribution follows as: 


R(t) 


u(r) = ROL) 9 


14 


From the momentum equation: 


Ou Ou 
= . +Uu dr 10 
p=|-p Cru) 





Using Equation 9 to evaluate the derivatives and integrating Equation 10 from the 


cavity wall, R, to infinity, we see that: 


RR+2R* pR? R* 


r 2 xr 
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P=p,+R 


Utilizing the boundary condition at the cavity wall where r = R and p(R) = 0, we see that 


Equation 11 becomes: 
O= p, + p-(RR +9, -R*) 12 


Solving Equation 12 using the initial condition of t = 0, R= 0, and R = Ro, we obtain 


the equation for the velocity increase of the collapsing bubble with radius: 


ee” R,) 
R’ P 0 Z 0 1 13 
30 R 
Integrating Equation 13 with respect to time yields the trajectory of the collapsing 
bubble: 
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The total collapse time, t, may be obtained from Equation 15 by integrating from zero to 


one. Thus, the total collapse time is: 


2 \/ 
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where B is the “Beta Function’, which is defined as: 


T(@)r 
B(a,B)= (a) (8) 16 
T'(a+ P) 
For the values of a, =5/6, 2 the collapse time becomes: 
hr 
r= aig 2 17 
Po 


Using Equation 17, it is possible to find the total time of collapse of a bubble of 
any given initial size. Realize, however, that this is only an approximate order of 
magnitude because compressibility effects become significant later on in the collapse and 
Equation 14 is no longer valid [Ref. 11]. The pressure distribution can be found from 


Equation 11 using Equation 12, which is expressed as: 


pa U abl (2 i }+(1-b) swere 6 =F 18 
Po 2 Po S a S RW) 


2. Converging Shock Waves 





Recalling the improved warhead design that was presented in Reference10, 
mathematical modeling can also be accomplished for converging shock waves. 
Previously, the collapse of bubbles in liquid was studied. Although the collapse of a 
bubble in a fluid medium does create a significant initial pressure wave, a greater initial 
pressure can be obtained by using a detonation device that focuses converging shock 
waves to a focal point that subsequently causes a rapid expansion. The asymptotic self 
similar solution was originally obtained by Guderley and later replicated by Butler. A 
detailed derivation of the solution to the problem of converging shock waves may be 


found in Reference 11 of this thesis. 
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HI. PREVIOUS EXPERIMENTS AND SIMULATIONS 
CONCERNING IMPLOSION PHENOMENA 


A. IMPLOSION EXPERIMENTS REGARDING STEEL CYLINDERS 

Long before copious computing power was readily available for educational and 
research oriented pursuits, experiments were completed in order to better understand the 
implosion phenomenon. During the early 1970s, a series of tests were completed by the 
Naval Ordnance Lab at White Oak, Maryland. During that time period, the implosion 
sounds of air-filled structures during collapse due to hydrostatic pressure had drawn the 
interest of the underwater acoustics community. For instance, implosions provided a safe 
and efficient way to study sound propagation in water [Ref. 13]. Also, studies conducted 
in this era by Reader and Chertock had shown that the pressure wave emitted by an 
imploding structure is negative while the interface is moving away from the water while 
the pressure wave is positive while it is accelerating toward the water. This principle is 
true whether the structure in question is glass walled, steel walled, or simply a cavity of 
air [Ref. 14]. 

1. Experimental Procedure 

The Naval Research Laboratory conducted implosion tests using steel cylindrical 
structures that were lowered or dropped into the ocean. The Naval Ordnance Laboratory 
recorded the acoustic signals produced as the steel tanks collapsed. The purpose of the 
experiment was to link the acoustic signals with the mode of collapse of the cylindrical 


structure. 


The test consisted of twenty-four tanks that contained air at atmospheric pressure. 
Each tank consisted of three compartments which were separated by rigid bulkheads. 
The center compartment was longer than the two outer compartments, which were of 
equal length. Thirteen tanks were controllably lowered into the ocean. After some type 
of collapse had been recorded, the tanks were raised back to the surface for inspection. 
For these tanks, external gages were mounted to record external pressures while an 


additional pressure gage was mounted internally at the center of the center compartment. 


As previously mentioned, some tanks were allowed to sink freely instead of being 


controllably lowered. Although this ruled out the possibility of recovery in most cases, 
17 


this type of testing allowed for more of the compartments to implode. One of these free 
sinking tanks had an internal pressure gage connected to a recorder by an expendable 
wire. In addition to the pressure gages, hydrophones were also used to collect acoustic 
signals emitted from the collapsing cylinders. Of note, the recording system performance 
was verified by recording signals from known explosive sources [Ref. 13]. 

2. Experimental Summary 

Inspection of implosion waveforms and the recovered tanks indicate that tanks of 
a similar material and geometry were damaged similarly and also created similar pressure 
waves. Generally speaking, the damage was uniform and the method of collapse was 
reproducible. An assumption was made that all cylinders collapsed because of 
hydrostatic loading by first forming an inward buckle. Further damage seemed to be a 
function of material properties, length to diameter ratio, wall thickness ratio, and the 


properties associated with the end cap structures. 


The steel tank implosion characteristics were compared against the glass bottle 
implosion characteristics which were studied in depth by Urick [Ref. 15]. In Urick’s 
experiments, the air-filled glass bottles propagated an initial negative pulse followed by 
diminishing pulsations. These later pulses were due to the gas bubble oscillation and are 
much like the gas bubble oscillations associated with underwater explosions. However, 
the initial pressure wave created by the implosion of the glass bottle is negative, while an 
initial pressure wave emitted from an explosion is positive. The following figure 
compares the pressure signals from one of Urick’s glass bottles to that emitted from an 


explosion. 
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Figure 8. Comparison of Explosion Pulse with Glass Bottle Implosion Pulse (From 
Ref.13) 


In comparison of the cylinder implosion signature with the glass bottle signatures, 
as well as comparisons against explosion wave forms, different waveforms were 
observed. The negative initial pulse of the glass bottle is due to implosion, while the 
oscillatory response is due to the gas bubble. The gas bubble induced response is much 
like the gas bubble oscillations in an explosion caused event. As for the steel cylinders, 
the negative initial pressure wave is attributed to the formation of an inward buckle. 
After the initial pressure wave, the response is largely dictated by the number and size of 
the buckles that form. The final collapse of the structure is indicated by a relatively large 
positive pressure pulse [Ref. 13]. Figure 9 displays the pressure time history at locations 


outside and inside of a steel cylinder. 
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Figure 9. Implosion of 26 Inch Steel Tank (From Ref. 13) 


B. SMALL-SCALE IMPLOSION TESTING OF GLASS AND ALUMINUM 
CYLINDERS 


Recently, implosion tests were conducted at Naval Undersea Warfare Center 
Division Newport. Four implosion experiments were conducted that consisted of 
cylindrical test samples made of aluminum (2) and borosilicate glass (2). The purpose of 
the experiment was to determine the effect of material properties on the rate of the 
implosion event. NAVSEA’s current method for assessing peak pressure created by an 
implosion event is a one dimensional, spherical bubble collapse model [Ref. 16]. 


NAVSEA’s model does not account for the structure that separates the air and water. 


The conclusions of the NUWC-NPT technical memo offer some insight into the 
events that occur during an implosion event. When an implosion occurs, the pressure 
wave is of the same magnitude at a constant radius from the event while varying the 
angle around the circumference of the cylinder. This effect was also true for an 
aluminum test sample, which does not collapse symmetrically. In addition, the rate and 


duration of the collapse is very much dependent on the material properties of the 
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structure. Also of interest, the experiments showed signs of irreversibility. Thus, losses 
or irreversibility is a characteristic of implosion phenomena [Ref. 16]. 
C; SIMULATIONS OF FLUID, GLASS, AND ALUMINUM MODELS 

Further research into the implosion phenomena has been completed at NUWC- 
NPT. In a more recent technical memo, personnel at NUWC-NPT document the 
development of a computational model that simulates the implosion of glass cylinders 
and spheres. Two models were constructed to analyze an implosion event. One model 
was a fluids-only model that simulated a low pressure gas surrounded by a high pressure 
water volume. The assumption was that at the initial time of t=0, the surrounding 
structure around the gas would uniformly disappear. The second model incorporated 
rigid body elements that at an initial time would begin to fail at a specified location and 
continue to fail at a specified “crack propagation rate.” For the second type of model, 
predictions were made for a sphere and four cylinders, each cylinder with a varying 


length to diameter ratio [Ref. 17]. 


Results from the simulations reveal variations due to length to diameter ratio and 
pressure probe location. Due to the distribution statement associated with Ref. 17, 


specifics can not be discussed in this thesis. 
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IV. STUDY OF IMPLOSION EVENT THROUGH MODELING 
AND SIMULATION 


A. INTRODUCTION TO DYSMAS 

Post World War II, the Navy completed many live fire UNDEX experiments on 
surface ships and submarines. Those tests provided useful data that helped create the 
basis for the empirical and semi-empirical equations used today in evaluation of peak 
pressures, decay rates, impulse, and other useful quantifications of UNDEX events. In 
addition to providing these useful tools, these live fire tests tested a vessel’s overall 


ability to withstand an attack. 


Even in modern times, similar experiments to those completed in the 1940s and 
1950s are done to test the combat survivability of modern war fighting ships. These tests 
are not only extremely expensive, but pose serious threats to the ship, the ship’s crew, 
and to the underwater environment. Recent tests included shock trials for the USS JOHN 
PAUL JONES and the USS WINSTON S. CHURCHILL and future shock trials are 
scheduled for the new LPD-17 class ship. 


Fortunately, a viable alternative to full scale shock trials is available. With 
modern advances in computing power, simulations of shock trials can be completed using 
the finite element method (FEM). At Naval Postgraduate School, data taken from the 
shock trials of both DDG-53 and DDG-81 was compared to simulations that were 
completed using LS-DYNA/USA, a finite element code. Close correlation to the actual 
shock trial was found. Further research into the simulation of these UNDEX events has 
provided strong evidence that finite element simulations are accurate approximations to 
ship response recorded during the ship shock trials. Figure 10 is a plot comparing the 


simulation results to the data taken from the live fire shock trials of DDG-81. 
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Figure 10. Midship Region Velocity Plot Comparison (DDG-81) 


Previously, LS-DYNA was used to simulate UNDEX events. However, a 
powerful Navy based alternative code, DYSMAS, has been released that offers 
advantages over the commercial LS-DYNA code. DYSMAS (Dynamic System 
Mechanics Advanced Simulation) is a coupled Eulerian-Lagrangian hydrocode that has 
been developed by Naval Surface Warfare Center, Indian Head Division. The code 
couples the GEMINI hydrocode with a version of DYNA-N. DYSMAS is especially 


useful for simulating the response of surface and subsurface structures to UNDEX events. 


GEMINI is a set of codes that is used to model the environment surrounding an 
UNDEX event. The environment can include the explosive, a liquid, or a gas that 
surrounds a ship, submarine, or other marine structure (however, DYSMAS is not limited 
to UNDEX events; it can also perform calculations for air blasts). The structural 
calculations are completed by a code such as DYNA-N. The Euler and structure codes 
are coupled using interface elements which exchange information after each 
computational step. Gemini supports both single processor and multiple processor 


computing [Ref. 12]. Figure 11 shows the relationship of the basic GEMINI components. 
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Figure 11. Basic Gemini Components (From Ref. 12) 


The process of using GEMINI is similar to most other FEM programs. Pre- 
processing and post-processing programs are needed before and after the calculation, 
respectively. In order to gain an understanding of the programs used in the GEMINI 
suite, a brief overview of the process is given below. A full input deck will be provided 
in Appendix C of this report. 

1. Pre-Processing Components of GEMINI 

In order to begin a calculation using DYSMAS, the environment must be 
established. This is done in the pre-processing portion of the GEMINI program. Setting 
up the surrounding environment for a calculation consists of setting up a coordinate grid, 
establishing an environment (water, air, explosive, etc.), and preparing for any parallel 


processing that is going to be done. 


The first step to setting up a simulation is to determine the coordinate system that 
will be used. GEMINI is capable of one, two, and three dimensional processing using 
spherical, cylindrical, or Cartesian coordinates. In addition the being able to use these 
three systems of coordinates, the field can later be rezoned in order to a higher dimension 


or different coordinate system. Once a coordinate system is chosen, the user can create a 
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unique flow field using GEMGRID. Using GEMGRID, the user can specify the size of 
the cells used in the calculation. Using GEMGRID is optional, if it is not used all cells in 
the calculation will have the same size. It allows the generation of the grid to be 
independent in each direction. GEMGRID is useful in cutting down simulation time, as 
cells away from the structure being studied can have an increased size (thus allowing for 


faster computation). 


Once GEMGRID has been run (if desired), the user then utilizes PREGEMINI in 
order to model the physical environment. In PREGEMINI, the user defines the water or 
air environment. In addition, the charge type, weight, and location are input in 
PREGEMINI. PREGEMINI is also used to ‘rezone’ previous flow fields (i.e., model a 
cylindrical 2-D calculation into a 3-D Cartesian calculation) or divide the computational 
domain into ‘subgrids’ that will allow for multiple processing. 

2. Computational Component of GEMINI 

Once the computational domain has been modeled using PREGEMINI, the main 
processing program can be used. GEMINI begins its calculation by referring back to the 
output files created by PREGEMINI. In addition to the PREGEMINI output files, the 
GEMINI program also reads its own user specified input deck. This user defined deck 
tells GEMINI is the simulation is a fluids only calculation or a coupled calculation 
(amongst other inputs). In addition to the type of calculation to be performed, other 
GEMINI inputs include cell, node, and element histories, termination parameters, restart 
parameters, and integration parameters. If a coupled calculation is being performed, the 
DYNA deck is also initialized with the execution of the GEMINI program. 

3. Post-Processing Components of GEMINI 

Once GEMINI (and DYNA if coupled calculation) has completed running, post- 
processing of the data can be performed. The GEMINI post-processing programs take 
the data that is output from GEMINI in binary form and translates that data into plot files 
that can be used in various visualization programs. There are two programs that are 
currently used to post-process the data: GEMHIS and GEMFIELD. GEMHIS is used to 
create a time history of points that were selected in the GEMINI input deck. These points 


can include structural nodes or elements or can be a point specified in the flow field. 
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There are a number of parameters that can be recorded, such as pressure, velocities, 


displacements, accelerations, density, forces, etc. 


In order to visualize the flow field, the GEMFIELD program is used. In 
GEMFIELD, the user specifies the subdomain of the simulation that is of interest. Then, 
GEMFIELD creates a time series of data records for the specified variables within that 


subdomain. GEMFIELD is useful when illustrations and animation are desired. 


Both GEMHIS and GEMFIELD are useful tools in visualizing the output data 
from GEMINI. The output from these two programs can currently be used in TecPlot 
and DysmasP applications. The output can also be put into plain ASCII if desired. In 
addition to specifying the variables to be plotted, GEMHIS and GEMFIELD also allow 
for the scaling of these variables. For example, this is helpful if the pressure output from 
GEMINI would want to be converted from dyne per square centimeter to kilopascals. 

B. SPHERICAL MODEL IMPLOSION SIMULATIONS 

As previously discussed, implosion models that have been previously studied 
include a fluids only model, a glass model, and a steel or aluminum cylinder model. The 
results from these tests indicate that a fluids only model will produce the greatest initial 
pressure wave. Intuitively, this is due to the uninterrupted inward flow of water creating 
what is similar to a converging shock wave. In the glass model, a crack propagation 
velocity and origin needed to be specified. In Steve Turner’s report, a velocity of 400m/s 
was chosen. Although this is a reasonable approximation, the velocity does vary and 
leads to the possibility of different magnitude shock waves being emitted from similar 
glass specimens. As shown in aforementioned experiments, the collapse of aluminum 
and steel cylinders are somewhat predictable, but also tend to give varying results 


between identical specimens in experimental tests. 


For the aforementioned reasons, a fluids only implosion modeling approach was 
utilized. Using a fluids only model, the pressure wave will be symmetric (in the case of 
an imploding spherical cavity) and will not depend on a crack propagation rate or 
location of initial damage. A fluids only model will also clearly depict bubble 


oscillations due to an implosion event. 
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In order to study the pressure wave due to an implosion, a spherical “air cavity” 
was created at a specified depth of water. At time t=0, the air cavity collapses upon itself 
due to the surrounding hydrostatic pressure. While the volume of the spherical cavity 
varied, the pressure probe locations were fixed relative to the radius of the cavity. Table 


1 indicates where the pressure probes were located in the simulation. 


Table 1. | Location of Pressure Probes with Respect to Cavity Radius 


Pressure Probe # X Distance Z Distance Radial Distance 
(m) (m) (m) 
0.0141 
|__| 1.414R 


3.536 R 
7.071 R 
14.142 R 





The probe locations are graphically illustrated in Figure 12. Again, notice that the 
locations of the probes are associated with constants based upon the air cavity initial 


radius. 





Figure 12. Pressure Probe Locations for Spherical Air Cavity Model 
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1. DYSMAS Simulations of Spherical Cavity Implosion 

As illustrated in Figure 13, the pressure decreased as the radial distance from the 
center of the implosion increased. Additionally, subsequent bubble oscillations create 
additional pressure pulses to a diminishing degree. This type of behavior is similar to 
that of an explosion event. The flow field due to a spherical implosion is illustrated in 


Appendix A. 
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Figure 13. Spherical Cavity Implosion of 2 Meter Radius Pressure History 


From the above figure, it seems apparent that the pressure wave created by the 
implosion acts very much like that of an explosion. The bubble pulses are similar and 
decrease as time progresses. Also, the pressure pulses that occur have less magnitude 
than the initial pressure wave. Also notice that as the distance from the center of the 
imploded volume is increased, the pressure wave magnitude is decreased. Indeed, from 
the above figure, it is apparent that the pressure wave created from an imploded spherical 
shape is influenced by the initial volume of the implosion and dissipates as the distance 


from the implosion is increased. 


The following figures are pressure history plots of implosion events. The volume 
of the implosion is varied. Recall that the probe locations are normalized by the distance 
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from the implosion center divided by the initial radius of the imploded gaseous sphere. 


The figures are scaled such that the initial pressure wave is exemplified. 
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Figure 14. Initial Pressure Wave Due to Spherical Cavity Implosion of 0.5 Meter 
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Figure 15. Initial Pressure Wave Due to Spherical Cavity Implosion of 0.75 Meter 
Radius 
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Figure 16. Initial Pressure Wave Due to Spherical Cavity Implosion of 1.0 Meter 
Radius 
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Figure 17. Initial Pressure Wave Due to Spherical Cavity Implosion of 1.25 Meter 


Radius 
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Figure 18. Initial Pressure Wave Due to Spherical Cavity Implosion of 1.5 Meter 
Radius 
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Figure 19. Initial Pressure Wave Due to Spherical Cavity Implosion of 1.75 Meter 
Radius 
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Figure 20. Initial Pressure Wave Due to Spherical Cavity Implosion of 2.0 Meter 
Radius 


Figures 14-20 have a similar pattern. Notice that as the volume of the imploding 
cavity increased, the time to the initial peak increased. This is intuitive, reasoning that 
more time would lapse before the larger volume would collapse upon itself. In addition 
to the time to initial peak increasing, the magnitude of the initial pressure wave also 
increases as the volume is increased. Of interest, the pressure recordings at the probes 
maintain the same relative trend regardless of initial volume. For example, notice the 
pressure history of probes one and three. In each of the simulations, the pressure time 
history of probe one generally cuts through the initial pressure wave peak of pressure 
probe three. 

2. Relation of Peak Pressure at Probe Locations and Implosion Volume 

As shown in the above data plots, the implosion pressure wave is influenced by 
the initial volume of the implosion. In order to determine how heavily the volume of the 
implosion influences the pressure wave, the pressure histories were plotted by individual 
probe. The data was then given a trend line that best described the data. The following 
figures display those data points and associated trend lines. Notice the good correlation 


between the data and the fit curve, as indicated by the R” value. 
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Figure 21. Peak Pressure Wave at Probe 1 as Function of Volume: Spherical Cavity 
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Figure 22. Peak Pressure Wave at Probe 2 as Function of Volume: Spherical Cavity 
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Probe 3 Peak Pressure P = 4.3038\0'%3 
R’ = 0.9665 
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Figure 23. Peak Pressure Wave at Probe 3 as Function of Volume: Spherical Cavity 
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Figure 24. Peak Pressure Wave at Probe 4 as Function of Volume: Spherical Cavity 
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Probe 5 Peak Pressure P = 2.274055 
R’ = 0.9021 
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Figure 25. Peak Pressure Wave at Probe 5 as Function of Volume: Spherical Cavity 


There are similarities between the probe histories. Notice the equations that best 
fit the data. Each of the equations is a function of volume to a power, preceded by a 
coefficient. Two assumptions should be noted. Acknowledging the fact that the pressure 
at probe one may be regarded as high in this comparison due to the short radial distance 
from the center of the implosion and, on the same token, the pressure history at probe five 
may be construed as low due to its relative further distance from the center of the 
implosion, these points (and their respective data fit trend lines) should weigh less 
heavily in the general formulation of the impact of volume to pressure. Thus, the 


0.132 The constant ahead of the volume 


influence of the volume seems to be that of V 
term is a function of the radial distance of that probe from the center of the implosion. 

3. Relation of Peak Pressure and Distance from Center of Implosion 

As mentioned above, the radial distance from the center of the implosion also 
impacts the peak pressure wave. In Figure 26, the peak pressure is plotted against 
normalized distance while the volume of the imploding cylinder and depth are held 
constant. Cleary shown is a dependence on radial distance, defined as (R/r)?”, Again, in 


this case the varying coefficients are due to the initial volumes of the imploding spherical 


cavities. 
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Figure 26. Peak Pressure as Function of Normalized Radial Distance 
4. Relation of Peak Pressure and Depth of Implosion 


Intuitively, as depth is increased the pressure wave created by an implosion will 
increase. This is due to the pressure difference between the pressure of atmospheric air 
inside the spherical cavity and the hydrostatic pressure of water outside of the cavity. 
Thus, not only should the peak pressure created by an implosion increase with increasing 
depth, it should increase linearly. In order to test this hypothesis, simulations were 
conducted that placed an equal volume sphere at depths of 38.1, 76.2, 152.4, and 304.8 
meters. Figure 27 displays the results of these simulations. Indeed, the increase of peak 
pressure is linear with respect to depth. However, of note is the trend that as the distance 
from the center of the implosion is increased (denoted by increasing probe numbers) the 


less of an effect depth has on the magnitude of the peak pressure. 
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Figure 27. Effect of Depth at Which Implosion Occurs on Peak Pressure 
3 Determining an Approximation to the Pressure Wave Created by an 


Imploding Spherical Cavity 

DYSMAS is a powerful modeling tool that has the capability to accurately predict 
responses due to UNDEX events. The results from the preceding tests demonstrated 
relations between the cavity’s volume and depth, the distance from the center of the 
implosion, and the subsequent pressure wave that was created by the imploding spherical 
air cavity. Using these similarities and relations, equations were derived that predict the 
peak pressure at a point due to an imploding spherical air cavity. Of course, the peak 
pressure is a function of depth, volume, and radial distance from the center of the 
implosion. In addition to presenting these equations, the error between the approximation 


and DYSMAS provided data is provided. 


First, the peak pressure will be determined by holding the depth constant (this is 
due to the fact that the power of the implosion is driven by a difference is pressure, which 
is constant for our hydrostatic case) while showing the influence of distance and volume. 


The equation should be of the following form: 
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P., =CV"R" 19 


where in this case, C is a constant that will later be associated with the depth of the 
implosion. From the data collected from the DYSMAS simulations, it has been 
determined that the radial exponent, n, is a constant value of n = -0.43 at a depth of 152.4 
meters. From Figure 27, it is shown that as the radial distance from the center of the 
implosion is increased, the rate at which pressure changes with respect to depth is 
decreased. Thus, the exponent of the radial distance term, n, is not a constant and varies 


with respect to depth. The exponent n is of the form: 


pa lTd_, OL z 
11430 300 


where d is the depth in meters. 


In addition to the varying value of n with respect to depth, the exponent of the 
volume, m, is not constant. From the collected data, it seems that as the volume of the 
implosion is increased, the amount that the volume influences the magnitude of the peak 


pressure wave increases. Thus, the exponent of the volume, m, is of the form: 


for r<1.25 m=(0.23 21 
r>1.25 tase 22 
15 75 


Lastly, the constant value must be determined. The constant is a function of 


depth. However, it is not linearly dependent on depth. The constant, C, is of the form: 
C =8.33E °d* + 0.025674d + 0.087759 23 
where again d is depth in meters. 


With the terms being defined, the approximation of the peak pressure created by 


an imploding spherical air cavity is defined as: 
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P.,. =(8.33E °d’ +0.025674d +0.087759)- oan ee 24 
R vA 14307 00 


where Prax is in MPa. The depth and volume should be entered in meters and meters’, 


respectively. 


This approximating equation is based on the data provided from DYSMAS 
simulations. The approximating equation provides good correlation to the data outside of 
the initial radius of the imploding spherical cavity. For example, if the spherical air 
cavity has an original radius of 1.0 meter, the approximation provides good correlation 
(under 10% error) outside of 2.0 meters away from the center of the implosion. Near the 
initial radius of the imploding sphere, the approximation under predicts the peak pressure 
by about 27%. With this in mind, the appropriate adjustment can be made for 
calculations at this range. Although the approximating equation does not correlate well 
within the original radius of the sphere, the damage done by an implosion event will 
likely occur outside of that radius. For example, a parent submarine hull would not be 
within the radius of the spherical air cavity, but it would be located close by. The 
approximating equation would predict the pressure wave that would impact that parent 


hull with relative certainty. 


While the approximating equation correlates well outside of the original spherical 
radius, it also does an adequate job at predicting the peak pressure at the center of the 
imploding volume. Most error is seen at the original radius of the spherical cavity. 


Figure 28 displays the correlating data for spherical implosion of 1.0 meter radius. 
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Figure 28. DYSMAS Data and Approximate Correlation: Spherical Cavity of 1.0 


Meter Radius 


The following table displays the approximated peak pressure as well as the error 
compared to the DYSMAS data. Notice that the percent error is at the most near the 


initial radius of the spherical air cavity. 


Table 2. Approximated Values of Peak Pressure and Associated Percent Error: Spherical 
Cavity of 1.0 Meter Radius 


R (m) Approximated Pmax (MPa) Error 
0.014 52.06 
1.414 28.1% 


3.536 
7071 
14.142 





Additional figures displaying the correlation between the DYSMAS data and 


approximated values are shown on the following pages. 
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Figure 29. DYSMAS Data and Approximate Correlation: Spherical Cavity of 0.5 


Meter Radius 


Table 3. Approximated Values of Peak Pressure and Associated Percent Error: Spherical 
Cavity of 0.5 Meter Radius 


R (m) Approximated Pmax (MPa) Error 


0.014 32.27 12.1% 
o7o7 | 600 254% 


1.768 
3.536 
7.071 
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Figure 30. DYSMAS Data and Approximate Correlation: Spherical Cavity of 1.5 


Meter Radius 


Table 4. | Approximated Values of Peak Pressure and Associated Percent Error: Spherical 
Cavity of 1.5 Meter Radius 
R (m) Approximated Pmax (MPa) Error 


10.607 
21.213 





This approximate equation is based on data from DYSMAS simulations. It is 
limited to the peak pressure calculation of spherical air cavities (at atmospheric pressure) 
placed in a hydrostatic field (water) with no structure. As the volume of the imploding 
sphere increases, the implosion becomes more unpredictable. Thus, this approximation 
should be limited to spherical air cavities with an initial radius of 2.0 meters or less. 

C. IMPLOSION OF CYLINDRICAL FLUID MODEL 

In addition to testing spherical air filled cavities, cylindrical air cavities were also 
modeled. Comparisons between spherical and cylindrical cavities of identical volume 
were studied. The probe locations for the cylindrical air cavity models are described in 
the table below and visually depicted in the following figure. Notice that the length to 


diameter ratio of the cylinders is, L/D = 2. 
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Table 5. Location of Probes for Cylindrical Air Cavity Tests 





Figure 31. Probe Locations of Cylindrical Air Cavity Tests, Quarter Model Utilizing 
Symmetry 


Using DYSMAS, two scenarios were modeled and simulated. In the first 
scenario, an air filled cylindrical cavity with L/D = 2 was modeled with the equivalent 
volume of a sphere of radius 0.5 meters. The second scenario consisted of an air filled 
cylindrical cavity of L/D = 2 with the equivalent volume of a sphere of 0.75 meter radius. 
The results were compared to the air filled spherical cavity results that were previously 
obtained in order to determine the difference in peak pressure magnitudes. Similar 


volumes were compared at equal distances from the center of the implosion. 


First, a comparison study considering the first scenario will be presented. Recall 


that the volume of the imploding sphere and cylinder are the same (volume equivalent to 
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a sphere of 0.5 meter radius). Figure 32 shows the pressure history of Probe 1, the probe 


at the center of both the cylindrically and spherically shaped implosion cavities. 
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Figure 32. Probe | Pressure History Comparison of Spherical and Cylindrical 
Implosions 


From the pressure history, it is evident that the implosion of the cylinder takes a 
longer period of time than the implosion of the spherical cavity. In addition, the peak 
pressure at this distance for the cylinder is nearly twice that of the peak pressure caused 
by the spherical implosion. The period of the bubble oscillation is also much greater in 
the case of the cylindrical implosion. Below is a table that compares the peak pressures 
caused by a cylindrical and spherical implosion. These implosions have the same volume 


and occur at the same depth. 
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Table 6. | Comparison of Peak Pressure Caused by Cylindrical and Spherical Implosions 


Distance from Pmax Pmax Greater Pressure 
Implosion Center Cylinder Sphere Pressure Ratio 
(m) (MPa) (MPa) 


0.0141 75.90 32.31 Cylinder 





Figure 33 displays the data from the table above. It is apparent that the peak 
pressure caused by the cylindrical implosion does not decline as quickly as the peak 
pressure caused by the spherical implosion. The pressure wave caused by the cylindrical 
implosion is approximately three times the magnitude of the pressure wave caused by the 


spherical implosion. 
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Figure 33. Comparison of Peak Pressure Caused by Cylindrical and Spherical 


Implosions 


Consider the second scenario, where the volume of both the sphere and cylinder 
are equal to that of a sphere with radius 0.75 meter. Again, the implosions occur at the 


same depth. The following table provides the peak pressure data for comparison. 
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Table 7. | Comparison of Peak Pressure Caused by Cylindrical and Spherical Implosions 


Distance from Pmax Pmax Greater Pressure 
Implosion Center Cylinder Sphere Pressure Ratio 
(m) (MPa) (MPa) 





Figure 34 displays the data from the above table. Notice that in two cases, 
distances of 5.3 and 26.5 meters, that the peak pressure is slightly elevated. This is due to 
the location of the probes. These probes are located perpendicular to the center of the 
cylinder, whereas the other probes (excluding probe 1) are located at an angle from the 
center of the imploding cylinder. As the cylindrical cavity implodes, the pressure wave is 
greatest along its axis and dissipates as the angle from the axis increases (of course, this 
is not the case with a spherical implosion due to symmetry). This phenomenon did not 


occur with the lesser volume cylinder. 
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Figure 34. Comparison of Peak Pressure Caused by Cylindrical and Spherical 


Implosions 


As shown in the modeled scenarios, cylindrical shaped implosions create a larger 
magnitude pressure wave than their respective equal volume spherical implosion 
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counterparts. The shape of the pressure wave is also dictated by the shape of the 
implosion. Clearly, the spherical implosion forms a symmetric pressure wave, while the 
cylindrically shaped implosion creates a pressure wave that is greatest along its axis of 
symmetry and decreases in magnitude as the angle from that axis increases. 
D. IMPLOSION OF ALUMINUM CYLINDER 

i Test Series 

Previously discussed was the data collected from fluids only models consisting of 
only water and air. Although those simulations offered insight into the behavior of an 
implosion event, a simulation is needed that would predict the pressure wave caused by 
the implosion of a structure. Experimental tests were completed in nearly forty years ago 
to determine the acoustical properties of implosion shock waves. More recently, 
simulations were conducted at NUWC-NPT to test the influence of length to diameter 
ratio of the cylinder, material construction, and initiation location on the magnitude of the 
implosion pressure wave. In those simulations, an initiating device was utilized that 
would create an initial deformation in the cylinder that would cause the cylinder to buckle 


and cause a pressure wave due to an implosion. 


Of interest is, instead of using an initiator to begin the implosion process, placing 
a small explosive charge near the cylinder and letting the pressure wave from the charge 
initiate the implosion of the cylinder. For these tests, the model will remain the same, 
allowing for a constant volume and identical material properties. However, the depth at 
which the cylinder is located (thus influencing the hydrostatic pressure surrounding the 
cylinder) will be varied. In addition, pressure probes will be located at specific locations 
to capture the effect that the cylindrical model has on the formation of the implosion 
pressure wave. Also, the charge weight will be varied to determine if decreasing the 
magnitude of the initiating pressure wave subsequently influences the pressure wave due 
to the implosion. The distance from the center of the charge to the outer hull of the 


cylinder will be held constant at 15.24 cm. 


A test matrix is given below that describes the types of tests that will be 


conducted in the simulations. 
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Table 8. Matrix for Structural Cylinder Implosion Simulations 


Charge Weight Depth 
(meters) 





The depths were calculated based on a percentage of critical buckling pressure. 
Tests one, two, and three are conducted at 50%, 75%, and 95% of buckling depth. The 


critical buckling pressure was calculated using the following: 


P=—. ~ 25 








where ¢ is wall thickness, r is mean radius, E is Young’s Modulus, and o, is the yield 
stress [Ref.18]. 

2. Computational and Structural Modeling 

The structural model was created using DYSMAS-P, a pre- and post-processing 
tool that complements the DYSMAS suite. The model dimensions are listed in the table 
below. 


Table 9. | Dimensions of Aluminum Cylinder 


Length 91.44 
45.72 
Wall Thickness 


The structural model is made of 4224 shell elements. The cylinder’s material 
property model includes an elastic/plastic behavior with strain rate dependence. Each 
shell element has three integration points through its thickness. The material model for 
AL6061-T6 was provided by NSWC-IH. Failure criteria were input to allow the deletion 


of failed elements. Figure 35 displays the cylindrical structural model. The endcap is 
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shown in red, while the body of the cylinder is shown in blue (for viewing clarity). Both 
structures are made of the same material. Notice that two planes of symmetry are used; 


half the longitudinal length and half of the circumferential distance are modeled. 


A, 


x 





Figure 35. Cylinder Structural Model 


The cylinder is placed at the predetermined depth and an explosive charge is 
placed 15.24 cm above the upper hull of the cylinder. In Figure 37, the charge would be 
placed at x=0, y=0, z=15.24+r (considering the shown endcap as L/2). A hydrostatic 
flowfield is placed around the model and the internal pressure of the cylinder is set to 


atmospheric. 


In order to decrease the computational time needed to complete the simulation, 
the initial pressure wave propagation was captured using a 2-D fluid pre-calc using 
GEMINI. Trap conditions were specified to halt the calculation once pressure or the flow 
field velocity change by a certain amount. Once the pre-calc completed, a rezone was 
performed in order to allow calculation in three dimensions. Figure 36 displays the 
trapped fluid flow condition for the 50% crush depth simulation. Again, the charge is a 
50 gram TNT explosive. 


50 


P (MPa) 
-109010 


200 
190 
180 
170 
-109015 160 
199020 
= 
oO 
NN 


-109025 


-109030 


-109035 





5 10 15 20 25 30 
x(cm) 


Figure 36. Fluid Pre-Calc for the 50% Crush Depth Cylinder Model: 50 g Charge 
Once the trap condition is achieved, the problem domain is rezoned into a three 


dimensional domain and the structure is added. Figure 37 shows the problem now 


rezoned into the three dimensional domain. Notice that the internal pressure of the 


cylinder is at atmospheric pressure. 
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Figure 37. Three Dimensional Rezone: Completed in PREGEMINI 


Once the re-zone has been completed in 3-D using PREGEMINI, GEMINI and 
DYNA_3D are used to calculate the responses of the structure and surrounding fluid. 
These simulations were only carried out to 8 milliseconds to capture the initial shock 
wave due to the implosion. A flow field displaying the pressure contour plots and 
structural response are shown in Appendix B. The response in Appendix B is for test 
four indicated in the test matrix, Table 8. Appendix C is an annotated input deck to run 


this simulation. 


In order to capture the effects of the pressure wave, cells have been denoted to 
record pressure histories throughout the simulation. The locations of these “pressure 


probes” are denoted in the following table. 
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Table 10. | Location of Pressure Probes; Coordinate Origin Considered at x = L/2, y=z=0 


7.0 | 0 | 1+25 25cm above; center 1 
[20d 0d 50450 om above: center —SS—S 
[3 [0 | 0 | +100 |f00.cm abovercenter—S—~S 
p40 80 *d om thwarts center SSCS 
ps [ 0 | +50 | 0 _]0 om athwart: center———SS 
Pe [0 re 100 [0 _|100 cm athwartcenter——SOSCS—~—~—~S 
[7 [0 | 0 | (+25) _]25 om below; center SS 
[eo | 0 | _-(r-¥50) _|0 om below: center SSS 
rs [0 | 0 | {r+ 100) |100 cm below; center SCS 
[i008] om above: quarterlengh ————SSC—S 
[attra 0 50] 50 om above: quarterlength SS 
peta 0 100] 100 om above: quarter length SS 
psa 50 *+d om athwart: quarter length SSS 
paar 50 | 0 ]0 om athwart: quarter length | 
p54 100 0] 100 om athwart, quarter length 
D6 ar 25) fe om below: quarter length SS 
[7a 0450) ]50 om below: quarter length SSS 
[isa 0 fr 100) 100 om below: quarter length SS 
[9 0 | 0 | 0 ____ internal; center cylinder SSCS 
P20 | t/a | 0 | 0 internal quarter length 





On the following page, a figure depicting the probe locations in accordance with 
the locations listed in Table 10 is displayed. Take note of the symmetry condition in this 


view. 
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Figure 38. Visual Depiction of Pressure Probe Placement Around Aluminum 
Cylinder 
3. Pressure Wave Behavior Around Cylinder 
In order to understand the behavior of the implosion shock wave at locations 
around the cylinder, the scenario of a 50 gram TNT charge with the model at 95% crush 


depth. The pressure probes used can be referenced in Table 10. 


First, pressure probe data will be examined from below the center of the cylinder. 
Recall that the small charge is placed above the cylinder, thus the fluid below the cylinder 
should not be influenced by the explosive’s pressure wave. Probes 7, 8, and 9 are 


displayed in the Figure 39. 
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Figure 39. Pressure History Below Center of Cylinder: Probes 7, 8, and 9 


Notice at early time (about 0.5 msec) that the ambient hydrostatic pressure is 
disturbed and pressure suddenly drops. This is due to the imploding cylinder; the fluid is 
rushing inward to fill the lower pressure volume. Then, an implosion occurs, creating a 


pressure wave that is four times that of hydrostatic pressure. 


Next, probes 4, 5, and 6 will be examined. These probes will display the pressure 
wave as it moves laterally from the imploding cylinder. Figure 40 displays the pressure 


histories of these probes. 
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Figure 40. Pressure History Athwart Center of the Cylinder: Probes 4, 5, and 6 


Notice the marked decrease as the distance from the implosion is increased. 
Because of the charge location (above the cylinder), the brunt of the implosion pressure 
wave travels in directions above and below the cylinder, while quickly diminishing 


laterally from the cylinder. 


Figure 41 displays the pressure history 50 cm around the cylinder at three 
positions; above, below, and athwart. Note that the pressure wave athwart of the cylinder 
is less. In addition, the initial peak of probe two (at time less than | msec) is attributed to 
the explosive pressure wave. The pressure wave caused by the implosion of the cylinder 


is also greatest in the direction toward the explosive charge (above the cylinder). 
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Figure 41. Pressure History Around the Center of Cylinder at Various Angles: 95% 


Crush Depth: Probes 2, 5, and 8 


Figure 42 displays the pressure history of probes also located 50cm from the 


cylinder at three positions, but the probes are located at a quarter of the cylinder’s length 


(not the center). 
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Figure 42. Pressure History 50cm from Cylinder at Quarter Length: 95% Crush 


Depth: Probes 11, 14, 17 


The difference in longitudinal position is not as influential as the angle from the 
cylinder. As shown in Figures 41 and 42 the pressure wave caused by the implosion is 
smaller at a position lateral to the cylinder. In addition, as the distance from the center of 
the cylinder increased, the magnitude of the pressure wave caused by the implosion also 


decreased slightly. 


Of interest is the pressure history for points that were initially inside of the 
cylinder. Probes 19 and 20 recorded the pressure history of this scenario. Probe 19 is 
located in the center of the cylinder while probe 20 is at the quarter length of the cylinder. 


Pressure histories for these two points are shown in Figure 43. 
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Figure 43. Pressure History of Probe Inside of Cylinder: 95% Crush Depth 


Notice that the pressure peaks at different times for the probes separated by 
22.86cm longitudinally. In addition, the peak pressure at the center of the cylinder has a 
slightly lower magnitude than the peak pressure at the quarter length of the cylinder. 

4. Comparison of Implosion Pressure Wave at Various Depths 

In order to determine the influence of depth on the implosion pressure wave, three 
simulations were conducted that varied the depth of the cylinder. The cylinder was 
placed in a hydrostatic field at 50%, 75% and 95% of hull crush depth. The charge size 
remained constant (a 50 gram TNT charge). Pressure histories were again collected at 


pressure probe locations, again as indicated in Table 12. 


Figure 44 displays the pressure wave 25 cm below the cylinder at various depths. 
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Figure 44. Probe 8 Pressure History at Various Depths 


Notice that as the depth is increased, the time until the implosion is accomplished 
is decreased. This is due to the increase of hydrostatic pressure outside of the cylinder. 
Of interest in this case is that the magnitude of the implosion pressure wave at 75% crush 


depth is greater than the 95% crush depth case. 
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Figure 45. Probe 19 Pressure History at Various Depths 


Figure 45 displays the pressure history at probe 19, which is positioned inside of 
the cylinder. Again, as the depth of the cylinder is increased the time elapsed until the 
implosion occurs is decreased. At this pressure probe, the magnitude of the implosion 


pressure wave increases as the depth is increased. 


Figure 46 displays the pressure history at Probe 3 at various depths. At this 
location, the time elapsed until the implosion occurs is decreased as the depth of the 
cylinder is increased. However, the magnitude of the pressure wave does not increase at 


the depth is increased. 


61 


607 


50 





40 





30 





20 





Pressure (MPa) 


10 























Oo 1 2 3 4 5 


Time (msec) 











50% Crush Depth 75% Crush Depth 95% Crush Depth 





Figure 46. Probe 3 Pressure History at Various Depths 


These simulations suggest that as the depth of the cylinder is increased, the time 
at which the cylinder implodes is decreased. There does appear to be a correlation 
between the depth of the cylinder and the magnitude of the pressure wave due to the 
implosion of the cylinder for the probes located at the center of the cylinder and above 
the cylinder. It is noted, however, that the magnitude and direction of the pressure waves 
are influenced by the implosion of the structure, which is not a process that occurs 
identically each time. 

5. Variance of Charge Mass and Influence on Implosion of Cylinder 

The mass of the charge influences the pressure wave emitted from an explosion 
event. Thus, by changing the mass of the charge, the pressure wave that will initiate the 
implosion of the cylinder will be affected. In order to study the effect of the mass of the 
explosive charge, three scenarios were modeled that placed a 12.5, 25, and 50 gram 
charge of TNT 15.24 cm above the center the cylinder. The pressure histories were again 


recorded using the probes located in accordance with Table 10. All the simulations were 
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conducted at 50% crush depth in order to mitigate any effects that depth may have on the 


implosion of the cylinder. 
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Figure 47. Probe 7 Pressure History: Effect of Varied Charge Mass 


Figure 47 displays the pressure history 25 cm below the center of the cylinder. At 
this location, the size of the charge does not seem to influence the time at which the 


implosion of the cylinder occurs. The effect of the charge size is inconclusive from this 


location’s data. 


Figure 48 displays the pressure history at the center of the cylinder. Again, there 
is no difference of time of implosion due to a variance in charge size. However, notice 


that the 50 gram charge does produce a lower implosion pressure wave at the center of 


the cylinder. 
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Figure 48. Probe 19 Pressure History: Effect of Varied Charge Mass 


Figure 49 displays the pressure history one meter above the cylinder. The initial 


peak is due to the explosive’s pressure wave (at about 0.5 msec). 
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Figure 49. Probe 3 Pressure History: Effect of Varied Charge Mass 
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From this set of three simulations, it is apparent that the charge size does not 
influence how quickly the implosion will occur. In these figures, time zero (t = 0) is the 
instant that the pressure wave has nearly reached the cylinder. From this point, the 3D 
simulation begins: the pressure wave interacts with the structure and the implosion 
occurs. The charge mass does not seem to have some impact on the magnitude of the 
pressure wave created by the implosion. Above the cylinder, the larger charge stimulates 
a larger implosion pressure wave response. However, that is not the case at the center 
and below the cylinder. Clearly shown in the data figures in this section, the larger 
charge weight did not necessarily produce a larger implosion pressure wave at all 
locations. However, in order to mitigate the influence of the explosive pressure wave, the 
smallest possible charge should be used to instigate the implosion of the cylinder. 

6. Effect of Charge Placement on the Implosion of a Cylinder 

In the previously modeled scenarios, the placement of the charge has been 
15.24cm above the cylinder. In the following modeled scenario, a 12.5 gram TNT charge 
will be placed 15.24 cm below the cylinder. The depth at which this implosion occurs is 


the same for the cases of the charge above and below the cylinder. 


The following figure displays the pressure histories of probes 100 cm above the 











cylinder. 
Explosion Above: Probe 3 (100cm Above Cylinder} 
Explosion Below Probe 3 (100cm Above Cylinder) 
2 3 
t (msec) 
Figure 50. Charge Placement Comparison: Pressure History of Probe 3 for Charge 


Placed Below and Above Charge 
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From Figure 50, it seems that the placement of the charge did not greatly alter the 


magnitude of the peak pressure wave. In addition, the general trend of the pressure 


histories is similar for both cases. The following figure displays the pressure history of 


probe 8, which is located 50cm below the cylinder. 
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Figure 51. Charge Placement Comparison: Pressure History of Probe 8 


Figure 51 shows that the peak pressure magnitude is not greatly influence by the 
placement of the charge. Again, as shown in both Figures 50 and 51, the trend of the 
pressure history is comparable between the two charge placements. 

The following figure shows the pressure history of probe 4, which is located 25cm 
athwart of the cylinder. Notice that the magnitudes of the pressure waves are similar and 


that the trend, as in the previous two figures, is similar between the cases with varied 


charge placement. 


66 





Charge Above: Probe 4 (25cm Athwart) 
Charge Below Probe 4 (25cm Athwart) 





P (MPa) 





2 
t (msec) 


Figure 52. Charge Placement Comparison: Pressure History of Probe 4 


The placement of the charge does not seem to have an effect on the pressure wave 
emitted from the imploding cylinder. As shown in Figures 50-52, the magnitude and 
general trend of the pressure wave are not impacted by the placement of the charge, 


whether the charge is above or below the cylinder. 


E. COMPARISON OF PRESSURE WAVES CREATED BY EXLOSION AND 
IMPLOSION 


Thus far, the pressure wave created by an implosion of a spherical cavity has been 
studied and an approximation method for calculating the peak pressure caused by the 
implosion has been accomplished. Of interest is the comparison of the pressure wave 
created by the implosion and the pressure wave created by a typical underwater 
explosion. Previously, the background and physics of UNDEX phenomena has been 
explained. Recall that the peak pressure at a distance, R, from a charge of a given type 


(TNT, HBX, etc) and weight can be calculated by using Equation 2: 


%" 
P =«,{ 2 





ma. 
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where K; and A; are properties of the charge type, W is the charge weight (Ibs), and R is 
the distance from the center of the charge (ft). Pinay is then determined (psi). Referring 
back to Figures 6 and 7, it is shown that the pressure wave created from an implosion can 
be greater than that of a similar explosion. This is accomplished by directing multiple 
detonators toward a central focal point (similar to a fluids only model in which the fluid 


begins moving toward the center of the cavity at the same rate from all directions). 


The following figure compares the peak pressure wave created by an explosive 
charge, a spherical fluids model implosion, and an imploding aluminum cylinder. The 
explosive is a 1.36 kg TNT charge. The spherical implosion pressure wave is based upon 
the approximation method discussed earlier. The parameters for the spherical implosion 
are a depth of 2073 meters and a diameter of 5 cm. The imploding cylinder occurs at 


2073 meters as well and the pressure readings are taken above the cylinder. 
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Figure 53. Comparison of Peak Pressure Wave Caused by Various UNDEX Events 


From Figure 53, it is apparent that the peak pressure wave created by an explosion 
and two modeled implosion scenarios are similar. Of note, however, is the relative small 
size of the charge and the extreme depth at which the implosions occur. The pressure 


wave from the implosion is a function of pressure difference (pressure outside of 
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imploding volume vs. the pressure inside the imploding volume). The pressure 
difference needs to be immense in order to develop the same peak pressure wave created 


by the small explosive charge. 
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V. FINAL REMARKS 


A. CONCLUSIONS 

It has been demonstrated through modeling and simulation that an implosion 
event does cause a peak pressure wave followed by subsequent pressure pulses similar to 
that of an underwater explosion. Previous research has proven that collapsing waves 
impinging on a central focal point can create a pressure wave that is greater than that of a 
typical explosion. This shaped charge weapon displayed in Figure 6 is an example of this 


principle. 


The shaped charge weapon initiates collapsing detonation waves that have a large 
pressure and velocity associated with them. In the models in this thesis, the driving force 
behind the collapsing wave was hydrostatic pressure. The difference between the internal 
pressure of the structure and the external hydrostatic pressure caused the implosion 
pressure wave. As shown through modeling and simulation, the magnitude of the 
pressure wave caused by an implosion of a small pressure vessel (with the driving force 
behind the initial collapsing pressure wave being hydrostatic pressure) creates a pressure 
wave that is comparable to that of a small underwater explosion. However, it seems that 
multiple detonators directing pressure waves inward toward a converging point are 
needed if an “implosion” pressure wave is to have a greater magnitude than a typical 


underwater explosive pressure wave. 


During experiments conducted by Naval Ordnance Laboratory, the collapse and 
implosion of cylinders under hydrostatic loading is shown to be somewhat repeatable. 
The depth at which implosion occurs is generally the same, but the shape of the imploded 
cylinder and the pressure pulse emitted from the implosion vary between occurrences. In 
the simulations conducted in this thesis, a cylinder was placed at a depth that was less 
than the crush depth of the cylinder. Then, a small charge was placed nearby the cylinder 
in order to initiate the implosion of the cylinder. Previous simulations have been 
conducted at Naval Undersea Warfare Center-Newport that initiated the implosion by 
placing a cylinder at subcritical hydrostatic pressure and initiating an implosion by giving 


the cylinder an initial deformation at a point. The pressure wave emitted from a “tapped” 
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cylinder implosion was symmetric about the entire model. However, when a small 
charge is used to initiate the implosion event, the pressure wave from the small charge 
interferes with the pressure wave caused by the subsequent implosion. This is shown by 
a non-symmetric pressure wave being emitted from the imploding cylinder as well as the 
general shape of the imploding cylinder. While a tapping of the cylinder initiates a 
controlled buckling process, the small charge creates buckling that collapses the entire 
cylinder upon itself. Also note, because of the interference caused by the explosive 
pressure wave, it is suggested that the smallest possible charge size be used in order to 
mitigate the influence of the explosive’s pressure wave. However, the placement of that 
small charge (whether below or above) does not have an effect on the pressure wave 
emitted from the implosion event. The use of a small charge to initiate the implosion 
process is recommended as a more realistic approach to modeling pressure vessel 
implosions. 
B. FURTHER STUDIES 

There is much work yet to be accomplished in the study of implosion phenomena. 
While this thesis has touched on the implosion of an unstiffened cylinder due to a small 
charge placed in close proximity to the cylinder, further work should be done on the 
implosion of stiffened cylinders. The implosion of a stiffened cylinder will be more 
complex, especially if the implosion is due to a placement of a small charge. The charge 
type and weight will be dependent upon the structure of the cylinder, and thus will be 


different for each structure studied. 


In addition, further studies should include simulations of shaped charge explosive 
devices. As documented in this thesis, a solid explosive charge outfitted with plural 
detonators that send multiple explosive pressure waves converging to a focus point can 
provide a marked increase in the magnitude of the pressure wave that is emitted from the 
weapon. The simulation of such an event should be accomplished. The future of mine 
warfare may include this type of warhead, and the effects of such a weapon should be 


studied. 
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APPENDIX A: 


SPHERICAL IMPLOSION FLOW FIELD 
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Figure 54. Spherical Implosion Flow Field: Frames 1-6 
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Figure 55. Spherical Implosion Flow Field: Frames 7-10 
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Figure 56. 


Spherical Implosion Flow Field: Frames 11-13 
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APPENDIX B: IMPLOSION OF ALUMINUM CYLINDER 
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Figure 57. Implosion of Aluminum Cylinder: Frames 1-6 
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Figure 58. Implosion of Aluminum Cylinder: Frames 7-12 
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Figure 59. 
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Implosion of Aluminum Cylinder: Frames 13-18 
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Figure 60. 
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Implosion of Aluminum Cylinder: Frames 19-24 
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Figure 61. 
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Implosion of Aluminum Cylinder: Frames 25-30 
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Figure 62. 
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Implosion of Aluminum Cylinder: Frames 31-36 
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APPENDIX C: GEMINI SUITE INPUT DECKS 


The following input deck comes from the implosion of an aluminum cylinder at 
95% crush depth using a 50 gram TNT charge. The following figure displays the flow 


path of this simulation. 








GEMGRID 2-D Cylindrical 
Coordinates 

















PREGEMINI |——~ | GEMINI 










































































| REZONE 
GEMGRID |——) PREGEMINI | ———- GEMINI 
+ 
DYNA 
3-D Cartesian Coordinates 
GEMFIELD GEMHIS 




















Post Processing 














Figure 63. Program Flow Path for Cylinder Implosion 


A. 2-DIMENSIONAL COORDINATE SYSTEM 
1. GEMGRID 























In this grid, the cells will be uniform in the X 
direction. The cells will be 0.25 cm in length. 
There will be 120 cells in this direction. 
<X BLOCK> 
+++++ 1,X mesh +++++ 
1 idatum 
om xdatum 
number xs xe ratio width 
NA 0.25 NA 1. 30 


<END X BLOCK> 
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#The Y BLOCK does not need to be input (only 2-D case). 









































<Y BLOCK> 
# ++++++ y mesh ++++++ 
fl; idatum 
O. ydatum 
#number xs xe ratio width 
1 Tey NA ‘ln: NA 
<END Y BLOCK> 
In this grid, the cells will be uniform in the Z 
direction. The cells will be 0.25 cm in height. 
There will be 128 cells in this direction. 
<Z BLOCK> 
t+++4 Z mesh +++++++ 
1 idatum 
-—207356 zdatum 
number xs xe ratio width 
NA 0.25 0.25 NA 32 
<END Z BLOCK> 
<POST> 
Standard # Half, Standard or double 


<END POST> 





2 PREGEMINI 





































































































<OPTIONS> 
START 
None # Modifications 
QO. # initial time (<0 then use value in restart file) 
-—980.665 # gravity constant(e.g. -98lcm/sec**2 ; z-direction) 
<END OPTIONS> 
<REZONE> section necessary for rezone runs only 
0 cell subdivision in each direction 
0 subdivision refine in mixed cells 
0. new flow field conservation scaling factor 
Os 12 0: location of old grid origin in new grid 
<END REZONE> 
<GRID> 
CYLINDRICAL # coordinates: CARTESIAN, CYLINDRICAL or SPHERICAL 
# cells cell width datum 
Xcells=120 dx=0. Xdatum=0. input number of cells from Grid.out 
Zcells=128 dz=0. dZdatum=0. input number of cells from Grid.out 
gemgrid/grid.asc calls the grid specified in GEMGRID 











<END GRID> 





<BOUNDARY CONDITIONS> 








# lower upper in each direction 
Xmin=wall Xmax=free b/c use of symmetry, the X min is fixed 
Zmin=wal] Zmax=free the 'free' allows for minimal reflection 


























<END BOUNDARY CONDITIONS> 
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<BLOCKED CELLS> 
imin imax jmin jmax kmin kmax 
<END BLOCKED CELLS> 
































































































































++++++ FLOW FIELD SPECIFICATION thts 
<EOS FILES> # Equation of State Files 
matl number matl file name ...... 
TNT D:\Seth\Tools\Eos_Lib\TNT.eos EOS file for 
TNT 
Tillwater D:\Seth\Tools\Eos_Lib\Tillwater.eos EOS file for 
Water 
Air D:\Seth\Tools\Eos_Lib\air.eos EOS file for 
Air 





<END EOS FILES> 





























<HYDROSTATIC FIELD> # creates a hydrostatic field 
zref=0. pref=1.e+6 zmax=max # Reference lin 
# matl name matl zmin matl int energy [EOS V1..Vn] 
air zmin=0. ei=eref 
Tillwater zmin=min ei=eref 





<END HYDROSTATIC FIELD> 




















<INITIAL STATES> 

#state matl g rho e p u Vv Ww 
1 TNT rho=rhoref ei=eref 

<END INITIAL STATES> 




















<FLOWFTIELD> 
option=hydrostatic # initializes hydrostatic field 
# option State 
# imin/max  jmin/max kmin/max state 
option=ball state=1 mass=50 x=0 z=-207341 # charge 
# parameters 
<END FLOWFIELD> 


























3. GEMINI 
<CASE> 
..\Pregemini\restart_000000.bin # start file (from PreGemini) 
IED Cylinder Pre-calc # TITLE (Limit 40 characters) 
Fluid # Fluid or Structure 
<END CASE> 










































































<TERMINATE> 

300 maximum step number 

0.003 termination time 

1.E-12 terminate if step size is less than this 
value 

# i/r <j j/y k/z variable change (+/-) This is a TRAP STOP 
x=30 j=l z=-207341 var=p dif=0.1 

x=30 j=l z=-207341 var=u dif=10. 

x=1 j=l z=-207355 var=p dif=0.1 

x=1 j=l z=-207355 var=w dif=10. 

x=1 j=l 2=-207327 var=p dif=0.1 

x=1 j=l 2=-207327 var=w dif=10. 

<END TERMINATE> 














85 


<INTEGRATION> 

-45 CFL safety factor 

05 initial step CFL factor 

2 limiter setting 

0 equalize after lagrange step (1l=yes, 0=no) 
1 equalize after remap step (l=yes, 0=no) 
1 

< 


























protect (O=off, 1 = on) 





END INTEGRATION> 








<CELL HISTORY> 
#X=7.62 zZ=-21.59 # 3.0 in r-dir 
<END CELL HISTORY> 























<CONTOUR PLOTS> 



















































































1000 step increments between plotting 
-5e-5 time increments between plotting 
<target plot steps> 
<end target plot steps> 
<target plot times> 
<end target plot times> 
<END CONTOUR PLOTS> 
<RESTART> 
100 restart tape write interval 
5 restart file save frequency 
<target restart steps> 
<end target restart steps> 
<target restart times> 
<end target restart times> 
<END RESTART> 
<TEXT OUTPUT> 
1000000 # step increments between printing [Gemini Only] 
1000. # time increments between printing [Gemini Only] 
standard # options: DEBUG or STANDARD, ELEMENTS is optional 

















# imin imax jmin jmax kmin kmax 
imin=l imax=1l jmin=1 jmax=1 kmin=10 kmax=1 
<END TEXT OUTPUT> 














KkKKKK RPND ***** 





B. 3-DIMENSIONAL COORDINATE SYSTEM 
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The 2-D portion of this problem consisted of setting up the environment and 
letting the charge explode. The 2-D portion was stopped when the pressure wave from 
the explosion was about to reach the cylinder structure. At this point, the calculation was 
‘rezoned’ into three dimensions in order to incorporate the structure. By doing a 2-D 
‘pre-calc’, the computational time was decreased. 


1. GEMGRID 


# This grid input is for a three dimensional Cartesian coordinate 
simulation. 





The X BLOCK displays the cell properties for the X direction. 
In this case, there are 4 lines of input. There are 20 cells 
of constant width (0.6 cm). Three are 35 cells that begin at a 
width of 0.65 cm and increase to 1.4 cm in width. There are 
5 cells that start with a width of 1.4 cm and end with 5.0 cm. 
There are 7 cells that begin with a width of 7.0 cm and end 


with a width of 62.8 cm. 
































<X BLOCK> 
++ ae ar mesh +++++ 
1 idatum 
0. xdatum 
number xs xe ratio width 
20 0.6 NA Ig NA 
35 0.65 1.4 NA NA 
2) 1.4 Dye.0) NA NA 
7 7.0 62.8 NA NA 
<END X BLOCK> 

















Like the X BLOCK, the Y BLOCK inputs indicate a nonuniform grid 
in the Y direction. 
<Y BLOCK> 
tee y mesh ++++++ 

idatum 

0. ydatum 
number xs xe ratio width 
23 O25 NA 1. NA 
12 1.0 NA ie NA 
5 1225 4.0 NA NA 
5 4.0 16.0 NA NA 
3 18.0 106.7 NA NA 


<END Y BLOCK> 





# Like the X BLOCK, the Z BLOCK inputs indicate a non-uniform grid 
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# in the Z direction. In addition, the zdatum is given to be at 


# Zmin. 


<Z BLOCK> 

# 

1 

-207469 

#number xs 
NA 
NA 
NA 
NA 
NA 
NA 


NOrROOFRN 
OO: OT Ci: © 


Standard 


+++4+4+4++ z mesh ++++4+++ 


idatum 
zdatum 
xe ratio width 
NA 1.0 50 
1.0 NA 15 
0.5 NA 50 
NA 1.0 15 
NA Lt 20, 50 
NA 1.0 70 


<END Z BLOCK> 





<POST> 
# Half, Standard or double 


<END POST> 





2. PREGEMINI 


<OPTIONS> 
REZONE 
BODYFILL 

0. 
-980.665 
<END OPTIONS> 























RE ZONE> 








04.505 
END REZONE> 


ANORUOUNUA 

















<SUBGRIDS> 

2 auto 

1 auto 

2 auto 

<END SUBGRIDS> 





<GRID> 
CARTESIAN # 
spherical) 
# cells 
Xcells=67 
Ycells=48 
Zcells=255 
grid\grid.asc 
<END GRID> 























:\Seth\Cylinder_Tests\ 





Rezone of the 2-D Precalc 





initial time (<0 then 


Allows a structure to b 


filled with property 
value in restart file) 

















gravity constant (e.g. 





























us 
=9 


8lem/sec**2 ; z-direction) 


section necessary for rezone runs only 
Explosion_95\Gemini2D\Gemini\restart_000080.bin 
cell subdivision in each direction 

subdivision refine in mixed cells 

new flow field conservation scaling factor 
location of old grid origin in new grid 








# Prepares for parallel processing 
# There are 2 subgrids created in X and Z direction 


coordinate (0 = Cartesian; 


cell width datum 


dx=0. Xdatum=0. 
dy=0. Ydatum=0. 
dz=0. Zdatum=0. 


# calls the Grid.out file 


<BOUNDARY CONDITIONS> 


# lower upper 


direction 
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= cylindrical; 2 = 


x or r cells 
y cells 
z cells 














































































































































































































Xmin=wall Xmax=free x or r 
Ymin=wall Ymax=free y 

Zmin=fr Zmax=fr Z 

<END BOUNDARY CONDITIONS> 

<BLOCKED CELLS> 

# imin imax jmin jmax kmin kmax 

<END BLOCKED CELLS> 

# ++4+4+++ FLOW FIELD SPECIFICATION +++++4+4 
<EOS FILES> 

# matl number matl file name ...... 
INT D:\Seth\Tools\Eos_Lib\INT.eos 
Tillwater D:\Seth\Tools\Eos_Lib\Tillwater.eos 
air D:\Seth\Tools\Eos_Lib\air.eos 
<END EOS FILES> 
<INITIAL STATES> 

#state matl g rho e u Vv Ww 
1 [NT rho=rhoref ei=eref 
2 Air rho=rhoref ei=eref 

S Tillwater rho=rhoref ei=eref 
<END INITIAL STATES> 
<BURN> 

# unburned to burned Det.time[s] x0 yO z0 
#Pentolitesolid to Pentolite Tbourn=0. x=0. y=0. z=-21.59 
<END BURN> 
<HYDROSTATIC FIELD> 

zref=0. pref=1.e+6 zmax=max # Reference lin 
# matl name matl zmin matl int energy [EOS V1..Vn] 

air zmin=0. ei=eref 
Tillwater zmin=min ei=eref 

<END HYDROSTATIC FIELD> 


<F 
# 
Op 


AOWF TELD> 
option State 
tion=hydrostati 











<E 








ND FLOWFIELD> 











Cc 
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The structure that is filled with a property is input in this 
section. 
The interfac lements and node file are needed. In addition, the 
location of the origin of the structure needs to be located in the 
flow field coordinate system. 


<BODY> 











Sk oe HR OE HE 



































Cylinder.ifa (15) interface element file (format identifer) 

Cylinder.nod node file 

0. O. —207326.2 #location of body coordinate system in 

fluid mesh 

1 Ox 0 vector along body coord sys x axis in fluid 
coord 

0 1. 0 vector along body coord sys y axis in fluid 
coord 

0 O. 1 vector along body coord sys z axis in fluid 
coord 

<END BODY> 





In BODYFILL, seeds need to be created that will fill the above 





mentioned 

structur These seeds will be the starting point of the fill, and 
will 

continue until a computational boundary, structural boundary, or 








subgrid boundary is reached. 


<BODYFILL> 

state=2 x=20. y=15.0 Z=-207380. 
state=2 x=40. y=15.0 Z=-207380. 
state=2 K=22.5 y=2.0 Z=-207357 « 
state=2 x=40. y=2.0 Z=-207357. 


<END BODYFILL> 





E 





EXT OUTPUT> 

# imin imax  jmin jmax 
imin=1l imax=1 jmin=1 
<END TEXT OUTPUT> 


kmin kmax 


jmax=1 kmin=1 kmax=1 




















a 


<CASE> 
PreGemini_Rezone\rezone_000000_xxx.bin 
Aluminum Cylinder Implosion 
characters) 

Coupled 

<END CASE> 


GEMINI 








file 
(Limit 40 

















or Coupled 











<TERMINAT 





Ez 
Vv 
































999999 maximum step number 

-008 termination time 

1.E-12 terminate if step size is less than this 
value 

# i/r 3/y k/z variable change(+/-) 

# 115 80 55 u 1.0% 

<END TERMINATE> 

<INTEGRATION> 
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45 




































































CFL safety factor 







































































































































































05 initial step CFL factor 
Zi. limiter setting 
0 equalize after lagrange step (l=yes, 0=no) 
1 equalize after remap step (l=yes, 0=no) 
ab protect (O=off, 1 = on) 
<END INTEGRATION> 
<TEXT OUTPUT> 
1000 # step increments between printing [Gemini Only] 
1000. # time increments between printing [Gemini Only] 
standard # options: DEBUG or STANDARD, ELEMENTS is optional 
# imin imax jmin jmax kmin kmax 
1 1 1 1 1 
<END TEXT OUTPUT> 
# --+--history----- 
<CELL HISTORY> 
# oi 5 k 
# xX y Z 
x=0.1 y=0.1 Z=-207331 25cm above center 
x=0.1 y=0.1 Z=-207356 50 cm above center 
x=0.1 y=0.1 Z=-207306 100 cm above center 
x=0.1 y=0.1 Z=-207428 25 cm below center 
x=0.1 y=0.1 Z=-207453 50 cm below center 
x=0.1 y=0.1 Z=-207503 100 cm below center 
X=22.85 y=0.1 Z=-207331 # 25cm above 3/4 length 
X=22.85 y=0.1 Z=-207356 # 50 cm above 3/4 length 
K=22.85 y=0.1 Z=-207306 # 100 cm above 3/4 length 
x=22.85 y=0.1 Z=-207428 25 cm below 3/4 length 
X=22.85 y=0.1 Z=-207453 50 cm below 3/4 length 
X=22.85 y=0.1 Z=-207503 100 cm below 3/4 length 
x=0.1 y=48.3 Z=-207380 25cm athwart at center 
x=0.1 y=73.3 Z=-207380 50 cm athwart at center 
x=0.1 y=123.3 Z=-207380 # 100 cm athwart at center 
K=22.85 y=48.3 Z=-207380 25 cm athwart at 3/4 length 
X=22.85 y=73.3 Z=-207380 50 cm athwart at 3/4 length 
K=22.85 y=123.3 Z=-207380 100 cm athwart at 3/4 length 
x=0.1 y=0.1 Z=-207380 center of cylinder 
K=22.85 y=0.1 Z=-207380 center of cylinder at 3/4 length 
x=0.1 y=7.9375 z=-566.9 4" freefield gage 
x=0. y=16.51 Z=-573.29 6.5" freefield gage 
x=0. y=46.99 Z=-573.29 18.5" freefield gage 
x=0. y=0. Z=-574.5 contact gage at cylinder crown, axial 
dist=0 (ctr) 
x=11.43 y=0.1 Z=-574.5 contact gage at cyl. crown, 4.5" 
axially from ctr 
X=22.86 y=0.1 Z=-574.5 contact gage at cyl. crown, 9.0" 
axially from ctr 
x=0. y=11.66 z=-585.0 # contact gage at 30 deg. down. 
x=0. y=23.6172 z=-605.5 contact gage at mid-draft, 90 deg. 
x=0.1 y=0.1 Z=-636 contact gage at cylinder keel 
x=0. y=0.1 Z=-4 15 one meter above top 
x=0. y=0.1 Z=-525 half meter above top 
x=0.1 y=0.1 z=-605.5 # center of imploded cylinder 
<END CELL HISTORY> 
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<CONTOUR PLOTS> 


















































20000 step increments between plotting 
100.e-6 # time increments between plotting 
<target plot steps> 

1 

<end target plot steps> 

<target plot times> 

<end target plot times> 

<END CONTOUR PLOTS> 

# +++++ RESTART FILES +++++ 

<RESTART> 

1000 restart tape write interval 

5 restart file save frequency 














<target restart steps> 
<end target restart steps> 


<target restart times> 
<end target restart times> 
<END RESTART> 





















































# t+++++++++++++ Required for Coupled Runs Only ++++++++++++4 
<COUPLED> 

l.et+6 back pressure 

.true. Load Failed SWI elements 

0. 0. -—207326.2 # location of Body coordinate system in Fluid Mesh 
Ly 0:8, Ose vector along Body coordinate system x axis in Fluid 
# coordinates 

One Lye “Os vector along Body coordinate system y axis in Fluid 
# coordinates 

Oe 0° Tx vector along Body coordinate system y axis in Fluid 
# coordinates 

<END COUPLED> 














<ELEMENT HISTORY> 









































4224 contact gage at cylinder crown, axial dist=0 (ctr) 
4213 contact gage at cyl. crown, 4.5" axially from ctr 
4203 contact gage at cyl. crown, 9.0" axially from ctr 
3804 contact gage at 30 deg. down. 

2922 contact gage at mid-draft, 90 deg. 

1578 contact gage at cylinder keel 





<END ELEMENT HISTORY> 

















<NODE HISTORY> 
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3 # crown node 

1122 # movin around mid-length circumference 
1154 

1186 








1570 
2 # mid-draft node 
4 # keel node 
# 3 # crown node 
292 # nodes along crown moving toward rear 
294 
296 
298 
300 
302 
304 
306 
308 
312 
316 
320 
324 
328 
1604 
<END NODE HISTORY> 


4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 











4. GEMHIS 



































<Case> 

all Process number (3 digits, "all", or "none") 
Subdirectory for output files (must previously exist) 
(don't use "/", use "." for current working dir) 

none Series: Output file identifier (1-10 characters or 
"none" ) 
(Series cannot end with "b" for DysmasP files) 

1 IVerbose: Screen output (0=no, 1l=some, 2=lots) 

-1.000 1.00 tBeg,tEnd: Time window for output (sec) 

0.000 1.00 tOffset, tScale: Offset and scaling for time 
t <== tScale* (toffset+t) 

ascii FileType ("ascii" or "binary") 

Tecplot FileFormat ("plain", "Tecplot", "DysmasP") 














<End Case> 


<Files> 
1 # IFileSplitGbl: Output file type (Global 


variable files) 
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1= single file 
2= one file for each variable 
3 # IFileSplitPt: Output file type (Point variable files) 
1= single file 
2= one file for each variable 
3= one file for each history point 
4= one file for each variable and history point 
<End Files> 
<Options> 
1 # ChargeScale: Scaling for specie (EnergyInt, Energy, Mass, 
Volume, Radius) 
-5.00e7 # pRefImp: Impulse intensity ref pressure (dynes/cm*2) 
# (make <0 to use value from lst record in file) 
0 # IRmOverlap (flag for removing oldest data if there is 


a time overlap) 





<End Options> 









































































































































He HE HE 
j 
oe 
O) 

5 
o 
3 
i) 
B 


se 
se 
te 
se 5 
se 
ste 
se 
ste 
se © 
se 
se 
se 
se © 





<Undocumented Options> 
#GageSymSize=8. Size (cm) of symbols in file gages.dat 
<End Undocumented Options> 














<Global Variables> 

# Key MatNum VarOffset VarScale #Note: MatNum=0 is for cell 
avg 

radius 1 0. Lg 

<End Global Variables> 











<Point Variables> 





Key MatNum VarOffset VarScale #Note: MatNum=0 is for cell 
avg 

u 0 O. 

ava 0 QO. 

Ww 0 0. 1. 

Pp 0 0. 0.0001 # converts to kPa 
r 0 0. 

i 1 O. 

Gj 2 0. 

e 0 0. 

e 1 Os 

e 2 O%, 

eint 0 O. 

eint QO. 

eint 2 0. 

ekin 0 0. 

if Oss 

£ 2 0. 
eos 1 0. 














<End Point Variables> 
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<Body Global Variables> 





# Key VarOffset VarScale 
force-x 0. alee 
force-y 0. dLfy 
force-z 0. ds 
force 0 Ls 





<End Body Global Variables> 








<Body Element Variables> 

# Key VarOffset VarScale 
p 0: ae 

<End Body Element Variables> 








<Body Node Variables> 

# Key VarOffset VarScale 
force-x O.. 
force-y 
force-z 
force 





DOCOCOWOOCOOC 0 


<End Body Node Variables> 

















































































































































































































































































































The following is a Key to the Global Variables 
HHEHE HH HEEHE EEE HEHE EEE HEHEHE SH HHEHE HEE HH HEHE HH 
Variable: allowed MatNum (0...NMat) 
energykin any 
energyint energy mass any 
vol radius >0 
volcav >0 
# tet HHEHE HEHE SH HHP HEE HEH Het HE Het HE 
The following is a Key to the Point Variables 
# tt HHEHH HEE SH HHEHE HHT HEH Het HE Het HE 
Variable allowed MatNum (0...NMat) 
all = all point variables NA 
u = velocity in r,x-dir 0 
v = velocity in y,theta-dir 0 
w = velocity in z-dir 0 
p = pressure 0 
r = density any 
e = total energy any 
int= internal energy any 
kin= kinetic energy 0 
f = volume fraction >0 
i = impulse intensity 0 
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eos#= EOS variable number NA 
HEHEHE HH HH HH HEEE HE HEH Ha EH EH EH EH EH HH HH HE HE EE 
The following is a Key to the Body Global Variables 
HEHEHE HH HH HH HEE HE HEH HHH EH HE HEH EH EH HH HH HE HE 
force-x = Total body force in r,x-dir 
force-y = Total body force in y,theta-dir 
force-z = Total body force in z-dir 
force = Total body force magnitude 
1 = Impulse 
HEHEHE H # HHH EH EH EH EH EE HH EE HE HE EE EE HEH EE EE EHH 
The following is a Key to the Body Element Variables 
HEHE HEH # Ht HE HEH aE ae Ea EH HH HE HE EE EE HE HEH EE EE HH 
p = Element pressur 
i = Impulse intensity 
Hitt # Ha HE HEH EH EH EH HH HH HE HE EE HE EE HE HE HE EH HH 
The following is a Key to the Body Node Variables 
Hitt # HEHE HEH EH EH EH HH HH HE EE HE EE EE HE EE HE EEE HH 
fForce-x = Node force in r,x-dir 
force-y = Node force in y,theta-dir 
force-z = Node force in z-dir 
force = Node force magnitude 
u Node velocity in r,x-dir 
v = Node velocity in y,theta-dir 
w = Node velocity in z-dir 
vel = Node velocity magnitude 
x = Node position (r,x-dir) 


# y= Node position (y,theta-dir) 
# z= Node position (z-dir) 


3: GEMFIELD 


<Case> 
all 


none 


1 
-1.000 1.00 
0.000 1.00 


ascii 
Tecplot 
<End Case> 





<File Output> 
3 


5000 


























Process number (3 digits, "all", or "none") 
Subdirectory for output files (must previously exist) 
(don't use "/", use "." for current working dir) 
Series: Output file identifier 

(Series cannot end with "b" for DysmasP files) 
IVerbose: Screen output (O=no, l=some, 2=lots) 
tBeg,tEnd: Time window for output (sec) 

tOffset, tScale: Offset and scaling for time 


























t <== tScale* (tOffset+t) 
FileType ("ascii" or "binary") 
FileFormat ("plain", "Tecplot", "DysmasP") 






































IFileSplit: Output file type 

1= single file 

2= one file for each variable 

3= one file for each time step 

4= one file for each variable and time step 
NTimeRecLimit: Max number of time records to generate 
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<End File Output> 


<Time Step Targets> 
#1815 6000 
<End Time Step Targets> 





<Plot Variables> 

















Key MatNum VarOffset VarScale #Note: MatNum=0 is for cell 
avg 

al 0 0. deg 

al includes the next nine variables 
u 0 Os 

Vv 0 0. 

Ww 0 0. ie 

Pp 0 0. 0.0001 # Convert to kPa 
ay 0 Ose 

r 1 QO. 

xr 2 0. 

e 0 0. 

e QO. 

e 2 0. 

eint 0 0. 

eint 0. 

eint 2 0. 

ekin 0 0. 

£ 1 QO. 

£ 2 QO. 

£ 3 QO. 

stat 0 Q.. 














<End Plot Variables> 


<Subdomain> 

xMin=0. xMax=0. iDelta=1 

yMin=-99999., yMax=999999. jDelta=1 
zMin=-9999999., zMax=999999, kDelta=1 

0 #NPtsRay 

O. O. 0. # xOffset, yOffset, zOffset (cm) 
1.e0 # xyzScale [xyz] <== 

[xyz] Scale* ([xyz]Offset+[xyz]) 
<End Subdomain> 




































































































































































Htttt Hitt tH Ht tt Ht Hitt Ht # 
#The next 2 sections are optional 
Hitt et Hitt H Hitt Ht Hitt Ht # 
<Options> 
EF Show parallel subdomain boundaries (default "F") 
F Show failed body elements (default "T") 
(Use "T" to preserve correct element numbering in Tecplot) 








A 


End Options> 





A 


Undocumented Options> 
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#addgeo=1l nodel=1l node2=69 

#addgeo=1l nodel=0 node2=2 

#ShowBlocked=T #logical flag to show values in cells 
blocked by body 

<End Undocumented Options> 
























































































































































































































































# # # HH EH HH EH EH EH EH EH HH HH HE EE EE HE Hitt 
The following is a Key to the Plottable Variables 
# # # HHH EH EH EH EH EH EH HH HH HE EE EE EE Hitt 
Variable allowed MatNum (0...NMat) 
all = all point variables NA 
u = velocity in r,x-dir 0 
v = velocity in y,theta-dir 0 
w = velocity in z-dir 0 
p = pressure 0 
r = density any 
e = total energy any 
int= internal energy any 
kin= kinetic energy 0 
f = volume fraction >0 
eos#= EOS variable number # NA [currently, only for P- 
alpha EOS] 
stat= status of cell (active/blocked) 0 
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